function varargout = PlotPsychFun_GAPS(varargin)
% x = PlotPsychFun_GAPS(schedule)
% x = PlotPsychFun_GAPS(schedule,trialstate)
% x = PlotPsychFun_GAPS(schedule,trialstate,CurTIDX)
% x = PlotPsychFun_GAPS(ax,...
% [x,pdata] = PlotPsychFun_GAPS(...
% 
% Plot psychometric function for gap detection paradigm
% 
% sensitivity is point on the psychometric function determined at 50%
% response probability.
% 
% DJS (c) 2010

if nargin == 1
    schedule = varargin{1};
    ax = gca;
    trialstate = schedule.trialstate;
    CurTIDX = length(trialstate);
        
elseif nargin == 2
    if isnumeric(varargin{1})
        ax = varargin{1};
        schedule   = varargin{2};
        trialstate = schedule.trialstate;
    else
        ax = gca;
        schedule   = varargin{1};
        trialstate = varargin{2};
    end
    CurTIDX = length(schedule.trialstate);
elseif nargin == 3
    if isnumeric(varargin{1})
        ax = varargin{1};
        schedule   = varargin{2};
        trialstate = varargin{3};
        CurTIDX = length(trialstate);
    else
        ax = gca;
        schedule   = varargin{1};
        trialstate = varargin{2};
        CurTIDX    = varargin{3};
    end

elseif nargin == 4
    ax = varargin{1};
    schedule   = varargin{2};
    trialstate = varargin{3};
    CurTIDX    = varargin{4};
end

trialstate = reshape(trialstate,length(trialstate),1);

r = trialstate;

x = strcmpi(schedule.writeparams,'dBSPL');
g = schedule.trials(1:CurTIDX,x);
gaps = unique(schedule.trials(:,x));

x = strcmpi(schedule.writeparams,'trial_type');
t = schedule.trials(1:CurTIDX,x);


for i = 1:length(gaps)
    %     n1(i) = sum(t == 0 & g == gaps(i) & ismember(r,[-2,1])); %#ok<AGROW>
    ms(i) = sum(t == 0 & g == gaps(i) & r == -2);          %#ok<AGROW>
    ht(i) = sum(t == 0 & g == gaps(i) & r == 1);           %#ok<AGROW>
    fat(i) = sum(g == gaps(i) & r == -1);                  %#ok<AGROW>
end

ms(isnan(ms))   = 0;
ht(isnan(ht))   = 0;
fat(isnan(fat)) = 0;

n1 = ht + ms;
n2 = sum(t == 1);

% ms  = ms  ./ n1 * 100;
ht  = ht  ./ n1 * 100;          ht(gaps <= 0)  = [];
fat = fat ./ (n1+fat) * 100;    fat(gaps <= 0) = [];

cr = sum(r == 2) ./ n2 * 100;
fac = sum(ismember(r,[-3,-4])) ./ n2 * 100;

g = gaps(gaps > 0);

ht(isnan(ht))   = 0;
fat(isnan(fat)) = 0;
cr(isnan(cr))   = 0;
fac(isnan(fac)) = 0;

if nargout > 0
    % fit line and find 50% response
    for i = 1:length(ht)-1
        if ht(i) < 50 && ht(i + 1) >= 50
            break;
        end
    end

    if ht(i+1) < 50
        sensitivity = nan;
    else
        [p,s,mu] = polyfit([ht(i) ht(i+1)],[g(i) g(i+1)],1);
        sensitivity = polyval(p,50,[],mu);
    end
end

try
    plot(ax,g,ht,'-ob',g,fat, ...
        ':r','LineWidth',2,'MarkerSize',10);
    
    hold(ax,'on');
    
    x = [min(g) max(g)];
    
%     plot(ax,x,[cr cr],'--b',x,[fac fac],':r','LineWidth',0.5);
    plot(ax,x,[fac fac],'--r','LineWidth',0.5);
    
%     plot(ax,[sensitivity sensitivity],[-5 50],'--k',[min(g) sensitivity],[50 50],'--k', ...
%         sensitivity,50,'*k','MarkerSize',10);
    
    hold(ax,'off');
    
    set(ax,'YLim',[-5 105]);
    set(ax,'XLim',[min(g) max(g)]);
    
%     x = sort([g; sensitivity]);
    set(ax,'XTick',g,'XTickLabel',num2str(g,'%3.1f'));
    set(ax,'XScale','linear');
    
    xlabel(ax,'Stim Level (dB SPL)');
    ylabel(ax,'% Total');
catch
    
%     fprintf('Problem plotting...\n')
end

if nargout == 0, return; end

varargout{1} = sensitivity;
varargout{2} = [g,ht'];
varargout{3} = fac;





